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Abstract 


A method for calculating orbital shadow terminator points is presented. The 
current method employs the use of an iterative process which is used for an 
accurate determination of shadow points. This calculation methodology is 
required, since orbital perturbation effects can introduce large errors when a 
spacecraft orbits a planet in a high altitude and/ or highly elliptical orbit. To 
compensate for the required iteration methodology, all reference frame change 
definitions and calculations are performed with quaternions. Quaternion algebra 
significantly reduces the computational time required for the accurate 
determination of shadow terminator points. 

Introduction 

In an effort to enhance the analysis capabilities of the Thermal Branch of the 
Structures and Mechanics Division at the Johnson Space Center, the authors are 
currently developing an analysis tool that will help in the design of a spacecraft 
mission attitude timeline based on the definition of thermal and power 
constraints. An important consideration of this task is the understanding of the 
duration of shadow passage for any given mission timeline, for this 
phenomenon has a significant effect on the on-orbit thermal environment 
experienced by a spacecraft. 

The problem of calculating the shadow times of a spacecraft orbiting a body has 
been studied in depth by various methods. 1-10 Assuming that the celestial bodies 
are spherical in shape, a planet's shadow consists of two distinct conical 
projections: the umbra and the penumbra (Figure 1). For the most part, however, 
the umbral shadow has been treated as a cylindrical projection of the Earth, 1-5 
since it significantly simplifies the calculation methodology. This assumption is 
fairly accurate for low altitude circular orbits but may lead to significant 
terminator point calculation errors for high altitude or highly elliptical orbits, in 
addition to ignoring the penumbral effects. Other authors 6-10 have treated the 
conical shadow projections. Peckman 6 treated only the umbral cone, while the 
others 7-10 have treated the effects of both shadow regions. However, these 
analyses have ignored orbit perturbation effects during propagation of the orbit. 
For thermal environment calculations, solar motion and perturbations give rise 

to the variation in the (3 angle, and ignoring these perturbations may lead to 

incorrect calculation of the environment, especially at (3 angle extremes where the 
spacecraft may be experiencing an uneclipsed orbit. The consideration of 
perturbation effects is required for an accurate calculation of shadow passage time. 
Additionally, only Dreher 8 has accounted for the effect of refraction in the 
shadow passage time. 


1 




Figure 1. Schematic representation of a planet's shadow regions. 


The current method is considered as an enhancement of the Long model. 4 
Through the use of an iterative procedure, the shadow passage times are 
calculated considering conical projections of umbral and penumbral shadows and 
perturbation effects. The current method does not consider refraction effects. 

Shadow Analysis Methodology 

The size and shape of the umbral and penumbral shadow regions are mainly 
functions of the planet size, the size of the sun, and the distance between the two 
celestial bodies. The refractive effects of the planet's atmosphere, although not 
considered in this analysis, could also affect the shadow geometries and therefore 
passage times. 7 ' 8 The umbra region is characterized by the total blockage of the 
solar energy component and the penumbra region by the partial blockage of the 
sun disk by the planet. In this region, the component of solar heating varies 
between a zero value at the umbra terminator to full solar at the penumbra 
terminator point. 7 

The calculation of shadow terminator points will be defined from the projection 
of the spacecraft onto the shadow cones, and the definition of the vector that 
points from the center of the umbral cone to the spacecraft at the projected site. 
The location of the spacecraft is calculated in the Mean of 1950 (M50) reference 
frame. The transformation calculations are performed using quaternion algebra, 
which significantly accelerates computation time. 

Definition of Shadow Cone Surfaces 

As in all of the shadow analyses considered, 6-10 this method assumes the celestial 
bodies to be spherical in shape, therefore, producing a conical projection of the 
shadow regions. (Note: the effects of a nonspherical body are accounted for only 
in the sense that they perturb the orbit.) This allows for the description of the 
umbral cone simply by considering the planet and sun diameters and the 
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separation between them. Then from geometry (Figure 2), we can calculate the 
umbral geometry from 


Xu 


D pV. 

®,-D p ) 


and 


(1) 
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=sin 
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(2) 


In the same way, the penumbral cone geometry can be determined (Figure 3) 
from 


X P = 


p A-. 

<D,+D„) 


and 


a p =sin 



(3) 


(4) 





Figure 2. Representation of the umbral cone geometry. 
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Figure 3. Representation of the penumbral cone geometry. 


Definition of Umbra and Penumbra Terminator Parameters 

The definition of the shadow terminator points is accomplished by locating the 
umbral and penumbral cones terminators at the projected spacecraft location. The 
location of the spacecraft in the M50 reference frame is represented by the 
r MS0 vector and is deferred to a later section. Figure 4 depicts the vector r s which 
defines this projection. The projection vector is obtained from the dot product of 
the vectors r M50 and s, that is 


F , = ( F M50*s)s. (5) 

Note that a shadow terminator may only be found when (F M50 » s)<0, as 
previously recognized by others. 4 With the definition of the r s vector, a second 
vector, 8, can be defined from 


^~ F M50 F s. (6) 

The 5 vector represents the distance between the center of the umbral cone and 
the spacecraft, at the projection point. Note that for the simplified assumption of 

a cylindrical umbral shadow projection, if the magnitude of the 5 vector is less 
than the radius of the planet, the spacecraft is considered to be in the planet's 
shadow. In the same way, the shadow terminator is found when the magnitude 

of the 8 vector is equal to the planet's radius. 4 " 5 Note that this analysis does not 
consider the case of a spacecraft orbiting a planet beyond the apex of the umbral 
cone. This assumption, however, is justified by the small subtended angles 
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associated with the umbral cone shadow geometry, which locates the umbral cone 
apex at great distances from the center of the planet. 



Figure 4. Representation of the r s and 5 vectors. 


In addition to the 5 vector, the determination of the r s vector allows for the 
location of the shadow terminator points at the specific projected location. From 
Figures 5 and 6, the distances k and £ are defined as 

^ = (Zu-|r s |)tana u (7) 

and 

^=(Zp+|r s |)tan« P - ©) 


The parameter % represents the distance between the center of the umbral cone 
and the umbral cone terminator, at the projected spacecraft location. In the same 
way, k represents the distance between the center of the umbral cone and the 
penumbra terminator, at the projected spacecraft location. 


As evidenced, a simple comparison between the magnitude of the 8 vector and k 

and £ defines the shadow terminator points. Specifically, the following 
comparisons may be drawn: 


a) Shadow terminators may only be encountered when (r M50 • s)<0. 

b) However, the spacecraft will still be in sunlight if |S|>k. 

c) The spacecraft is in penumbra if £ < |s|<k. 
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d) The spacecraft is in umbra if 5|<^. 

e) A penumbra terminator point is found when |s|=k. 

f) An umbra terminator point is found when |5|=£. 



Figure 5. Location of the penumbral cone terminator at the projected spacecraft location. 



Figure 6. Location of the umbral cone terminator at the projected spacecraft location. 

Since the methodology only deals with the magnitudes of the 8 vector and the k 
and £, parameters, the determination of an entry or exit terminator point requires 
additional consideration. If the analysis is performed by advancing in eccentric 
anomaly 11 until the orbit is completed, then the following observations can be 
made: 


a) If at the beginning of the analysis S|>k and |s|>£/ the spacecraft is initially in 
the sunlight. The first terminator encounter, if any, must be a penumbra entry 
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point. The second terminator encounter may be either an umbra entry or a 
penumbra exit. If a penumbra exit is found, the analysis for this orbit is 
completed. If an umbra entry point is found, the third encounter must be an 
umbra exit, followed by a penumbra exit. Then, the orbit analysis is completed. 

If tpin is the time of penumbra entry, tpex is the time of penumbra exit, t u in is 
the time of umbra entry, and finally tuex is the umbra exit time, the time of 
shadow passage is determined as 

Time in umbra = tuex * tuin (9) 

and 

Time in penumbra = tpex - tuex + tuin - tpin- (1°) 

b) If at the beginning of the analysis 6 <k and 8 <\, the spacecraft is initially in 

umbral shadow. The first terminator encounter must be an umbra exit point, 
followed by a penumbra exit. After a period of sunlight, the penumbra entry 
point is found, followed by the umbra entry. The finding of all terminator 
points completes the analysis of the orbit. 

If tper is the period of the orbit, the time of shadow passage can be calculated as 
Time in umbra = tuex + tper - tuin (11) 

and 

Time in penumbra = tpex - tuex + tuin - tpin- (12) 

Since the initial problem time is only an offset, it does not appear in the time 
equations. 

c) If at the beginning of the analysis 8 <k and 8|> %, the spacecraft is initially in 

the penumbral region. Then the first terminator encounter can either be an 
umbra entry or a penumbra exit. If the penumbra exit is found, the analysis is 
completed for the orbit. If instead, the umbra entry terminator is found, then it 
must be followed by umbra exit, penumbra exit, and finally by a penumbra 
entry. 

The time of shadow passage can be determined from 

Time in umbra = t U ex - tuin (13) 

and 
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Time in penumbra = tpex + tper - tpin - Time in umbra (14) 

Finally, the time in sunlight is simply calculated as the total orbit period less 
the time spent in shadow. 

Determination of the Spacecraft Location 

The definition of the relationship that exists between two coordinate systems 
sharing a common origin is the basis for the calculation of the r M50 vector. This 
definition may be obtained from any of three mathematical treatments: mainly 
the Euler angle, the direction cosine matrix, and the quaternion. Until the 
development of the digital computer, the Euler angle definition was widely used 
due to its geometrical simplicity and clear visualization. 12 The digital 
computation advancement of the mid-1960s marked the beginning of the use of 
the direction cosine matrix treatment since the computation methodologies were 
more suitable for computer programming, particularly when successive 
transformations of a body with respect to a fixed reference frame were defined. 

The use of the direction cosine matrix methodology is still common today. The 
third mathematical treatment is through definition of the quaternion. 

Quaternions 

An infrequently used mathematical treatment of body transformations is the 
quaternion. This treatment was first devised by Sir William Rowan Hamilton 13 
in 1843. This approach makes use of Euler's theorem which states that any real 
transformation of one coordinate system with respect to a fixed reference system 
(sharing a common origin) can be described through a single rotation called 
principal rotation about a single axis called principal axis. 

A quaternion is a compact representation of a principal rotation about the 
principal axis and can be represented 12 " 17 as an ordered quadruple of real numbers 


qAB-tq1.q2.q3.q4l 

r f . f . f . f, 
qAB =tcos-,e x sin-,e y sin-,e z sin-], 

or expressed in vector form 13 as an addition of a scalar and a vector 

= scalar + vector 

qAB=qi+q2i+q3j+q4k. 

The definition of the quaternion is subject to the normality condition 


(15) 


( 16 ) 
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q, +q2+q3+q4-l- 


(17) 


The treatment of quaternions is much like the direction cosine matrix in that a 
successive order of transformations results in a total transformation quaternion 
that is obtained through successive quaternion multiplication. The advantage of 
using quaternions results from the reduced computational load associated with 
the calculation of the total transformation quaternion. However, the 
computational load savings are only realized when a conversion of the total 
transformation quaternion to the equivalent total transformation matrix is not 
required and when the interchangeability property of the quaternion 
multiplication 12 is used whenever possible. As seen from the representation of 
the quaternion, the definition the quaternion requires four elements versus nine 
elements required to define a direction cosine matrix. Therefore, an added benefit 
is immediately realized from the reduced computer memory requirements and 
number of elements that need to be manipulated. It is therefore recognized that 
by using quaternion algebra, the calculation of transformations required for the 
determination of shadow passage are significantly improved. This allows for the 
use of an iterative process without severely impacting computation time. 

Quaternion Algebra 

The advantages of quaternion manipulation for guidance and control were 
recognized early in the development of the Space Shuttle Orbiter. Several 
internal NASA publications 1517 describing the Shuttle onboard software 
manipulation of quaternions were an important consideration in the 
development of the quaternion algebra defined in this section. Most importantly, 
the quaternion multiplication and the vector transformation through a 
quaternion will be described. 

The vectorial definition of the quaternion allows for the development of 
quaternion algebra in the classical sense. If two quaternions, Q and P, are defined 
as 


Q-fqi^/Cb'CU) 


P=(p 1 ,p 2 ,p 3 ,p 4 ) 


(18) 


then the fundamental definitions 17 are 

a) Equality: Q=P, when and only when 

qi=Pi>q 2 = P2/q3=P3,a il dq4=P4 (19) 

b) Addition: 
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Q + P = (q 1 + p„q 2 + p 2 ,q 3 + p 3 ,q 4 + P 4 ) 


( 20 ) 


b) Subtraction: 

Q-P = (q 1 -pi,q 2 -p 2 ,q 3 -p 3 ,q 4 -p 4 ) 

c) Multiplication by a scalar: 

aQ = (aq 1 ,aq 2/ aq 3 ,aq 4 ) 

d) The quaternion product: 

QP = (q, + q 2 i + q 3 j + q 4 kXpj + p 2 i + p 3 j + p 4 k). 


( 21 ) 


( 22 ) 


(23) 


If we express the scalar part as S and the vector part as V, the product may be 
written as 

QP=(S q +V q )(S p +V p ). (24) 

Manipulating this expression we obtain 

QP=S q S p +S q V p +V q S p +(V q xV p )-(V Q •V P ). (25) 

This expression has been shown in the literature 12 ' 13 ' 17 in matrix form as 


QP = 


qi 

q 2 

q 3 

l_q 4 


-q 2 

q 3 

-q 4 

q 3 


-q 3 -q 4 
q 4 -q 3 
qi q 2 
-q 2 qi 



>r 


p 2 


Ps 


J>4. 


(26) 


If the order of the quaternion is reversed such that 

PQ=S p S q +S p V q + V p S q +(V P x V q )-(V p *V q ), 
the resultant matrix form is 


(27) 


PQ = 


Pi 

P 2 

P 3 

Pi 


-P 2 

Pi 

P 4 

-P 3 


-P 3 

-P 4 

Pi 

P 2 


-P 4 

P 3 

-P 2 

Pi 



'q/ 


q 2 


q 3 


.q4. 


(28) 
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If the resultant matrix forms of equations 26 and 28 are simplified, both will give 


(QP)i =( p Q)i =Piqi -p 2 q 2 -p 3 q 3 -p 4 q 4 
(qp ) 2 =(pq) 2 =p,q 2 +p 2 qi +p 3 q 4 -p 4 q 3 
(QP ) 3 = (PQ ) 3 =Piq 3 -p 2 q 4 +P 3 qi +P 4 q 2 
(qp ) 4 =(pq) 4 = p,q 4 +p 2 q 3 -p 3 q 2 + p 4 q 3 


showing the interchangeability of the quaternion multiplication. 

The 4X4 matrices of equations 26 and 28 are defined 12 as the quaternion matrices. 
Notice that the only difference between quaternion matrices in equations 26 and 
28 is the transmuted nature of the minor matrix of the first element. If we 
describe a transformation from the A reference frame to a B frame, and then from 
the B to the C frame, we can express the quaternion multiplication in a matrix 
form as 


Qc<-A - [M] c< _ b Qb<-A- 


(30) 


or 


Qc.A - [M]b<-A Qc«-b 


(31) 


where [M]‘ B< _ A is the transmuted quaternion matrix of the A->B transformation. If 
an additional transformation from C->D is imposed, then the total transforma- 
tion can be expressed in matrix form as 

Qd<-a = [M] D( _ c [M] Ct _ B Q b *_ a (31) 

or 

Qd^a = [M] d ^ c [M]^ a Q c ^ b (32) 


and 


In general 


Qd<-A — ML a [m] d ^ c q c ^ b . 

D <r-C [M] 

D<-C [M]La , 


(33) 


(34) 


is a property that can be extended to the product of any number of quaternions. 
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Vector Transformations 


Given, without proof, 15 " 18 a vector is transformed from one coordinate system to 
another by 


~ Qd«-a v a Qd«-a/ (35) 

where Qo<_ A is the conjugate of Q D< _ A . The conjugate of a quaternion represented, 
for example, by equation 15 is defined 15 " 18 as 


(36) 

The manipulation of equation 35 is through normal quaternion multiplication, 
where v A is treated as a quaternion with a zero scalar element. 

Motion in the Orbit Plane 

The motion in the orbit plane, as shown in Figure 7, is obtained in part from the 
polar equation of the ellipse 


a(l-e 2 ) 

1+ecosv' 


(37) 


which can also be expressed in terms of the eccentric anomaly 11 as 

r=a(l-ecosE). (38) 

In Figure 7, the subscripts sc and p refer to the spacecraft and planet coordinate 
systems, respectively. 

The relationship between the eccentric anomaly and the true anomaly 11 is given 

by 



v= 



tan — E. 
2 


The time of flight is calculated from Kepler's equation 11 ' 14 


M=E-esinE= 



(t-T), 


where T is the time of passage through pericenter. 


(39) 


(40) 
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Calculation of the r MS0 Vector 

As mentioned previously, the calculation of the r M50 vector is performed in the 
M50 reference frame. The description of this vector is shown in Figure 8, as 
defined by the Keplerian elements. 

Note that the orbit perturbation parameters are modeled as time variations in 
right ascension of the ascending node, Q, and argument of pericenter, (0. 



Figure 7. Polar coordinate definition of in-plane orbit parameters. 
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Figure 8. Definition of spacecraft location in the M50 reference frame. 


The r M50 vector is obtained from 


^M50 — r ‘(r p ^sp)/ (41) 

where r p ^ sp is the unit vector that points from the center of the planet to the 
spacecraft, and r is obtained from equation 37 or 38. The unit vector r is 
obtained from a simple sequence of transformations, mainly 

Q?M50 = Q(a>+v) Q; Qfl (42) 


and 


"?MS0 — Qm50 Qm50' 

where ? y is the reference vector {1,0,0}. 


(43) 
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In equation 42, the rotation sequence is obtained by first, a z-axis rotation 


Qn = [cos y , 0. , 0. , - sin y ] , 


( 44 ) 


followed by an x-axis rotation 


. Qi=[cos^,-sin^ / 0. / 0.] / 


(45) 


and finally by another z-axis rotation, which by equation 20 can be represented as 


r (o) + v) _ . (00 + v), 

Q(«+ v) = [c° s — i — , 0. , 0. , - sin ^ ; ]. 


(46) 


It should be noted by the reader that with little or no variation in right ascension 
of the ascending node and orbit inclination, the calculation Q, Q n of equation 42 
may be performed only once for each orbit analysis, thereby saving a significant 
amount of computational time. 


In quaternion matrix form, equation 42 may also be expressed as 


Q rM »=[ML +v) MQa (47) 

which by relations 32 and 33 can also be expressed as 

( 48 ) 


Once again, with little or no variation in Q and i, the matrix manipulation 
[M]q [M]‘ may be performed only once per orbit analysis. 

With the flexibility that quaternion algebra offers, the terminator calculation 
method employed here is well suited for considering the variation in H and co 
due to the effect of the J 2 oblateness perturbation. The variation of these two 
parameters will have the most pronounced effect on the calculation of the 

terminator entry and exit times. In such a case, the and co parameters are 
determined by the following expressions: 


£2(t) =Q 0 +Dt 


(49) 


and 


o}(t)=co 0 +d)t / 


(50) 
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where £2 C is the initial right ascension of the ascending node, and co 0 is the initial 
argument of pericenter. In equations 49 and 50, £2 and (b are obtained^ from 



(51) 


and 


• 3 t 




[5cos 2 (i)-l] 


(52) 


where req is the equatorial radius of the planet/ moon, p is the parameter or semi- 
latus rectum, and J2 is the oblateness perturbation coefficient. 


Iterative Methodology 


As mentioned in previous sections, the calculation of the shadow terminator 
points is solved by an iterative procedure. The iterative method presented here, 
however, is well-suited to inclusion of the perturbations while the orbit 
progresses. Hence, no prior determination of whether or not the solar motion 

and the rates of co and £2 significantly affect the terminator locations need be 
made. A standard bisection method, as described for example, by Cheney, et al., 19 
was adopted for this application. The bisection method proves to be the most 
effective calculation methodology, since the shadow terminator variables 


A p =(|S|-K) and A u =(5 -\) will have opposite signs in the time interval 

containing the terminator, and the shadow function described by the time 
variation of the shadow terminator variables is continuous. 


The analysis of an orbit is performed by advancing in eccentric anomaly. An 
analysis interval [A a ,A b ] / represented by 


A^A^DWl (53) 

and 

A,. k =A,. b [E»(t t )l, («) 


(<p is either p or u for penumbra or umbra, and tb>ta) will contain a terminator 
point if A > 0 and \ rb < 0. To find this terminator point, an intermediate 
analysis point is created as represented by 
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\,c=\,c [E c (t c )] 


(55) 


where 


tc=|(t a +t b ). (56) 

Ideally, a shadow terminator point if found at time tc (note that tb>tc>ta), if 
A q) C = 0. However, this is seldom the case. Instead, a sequence of time interval 
assessments must be performed according to the observation that if 

A<mV <0 ' (57) 

then the shadow terminator must be located between times t a and tc- Similarly, if 

A^A^O (58) 

the shadow terminator is then in the time interval [tc,tbl- A new interval 
assessment is performed by selecting the time interval containing the shadow 
terminator. This sequence is repeated until A v c - A, a < error or A b -A < error, 
where the error is a satisfactorily small number. 

Sample Cases 

The mathematical methodology described in the previous sections was coded in 
standard FORTRAN in support of the development of the Thermal Constraint 
Attitude Design System (TCADS) analysis package. TCADS will help in the design 
of a spacecraft mission attitude timeline based on the definition of thermal and 
power constraints. 

Five examples were chosen based on simulations that would exercise the 
capabilities of the methodology to the full extent. All simulations were performed 
with orbits about Earth. However, the method is applicable to any celestial body 
for which the necessary parameters are known. Physical constants used in these 
analyses were obtained from the Jet Propulsion Laboratory. 20 The examples are 
named as follows and covered individually in the following sections: 

•High Inclination, Low-Earth Circular Orbit 
•Sun Synchronous Orbit 
•High Inclination Elliptical Orbit 
•Geostationary Orbit 

•Combined Released and Radiation Effects Satellite (CRRES) 10 Case Comparison 
•Illustration of Conical versus Cylindrical Shadow Assumption 
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In the first five examples, the shadow terminator analysis was performed for 2 
Earth years. The results are presented in plots of beta angle, p, variation versus 
time, and variation in percentage of orbit period spent in umbral and penumbral 
shadow as a function of time. The last example provides an assessment of time 
spent in Earth's shadow for both the cylindrical assumption 5 and the conical 
assumption using the current methodology. 

The beta angle, P, is the angle between solar vector, s, and its projection onto the 
orbit plane, and it is used in this analysis since it provides the thermal engineer a 
’ means of assessing the thermal environment experienced by an orbiting 

spacecraft. P is given by 


P=sin _1 (o*s), (59) 

where 6 is the orbit normal vector and s is the unit solar vector. These unit 
vectors may be expressed by 


f cos(D 


§= 


sin(T) cos(e) > 


sin(r) sin(e)J 


( 60 ) 


and 


sin(Q) 

■cos(Q) sin(i) 
cos(i) 


( 61 ) 


where T is the ecliptic solar longitude and e is the obliquity of the ecliptic. As a 
planet moves about the sun, T will vary from 0 to 2ic. Additionally, the 

perturbation in Q, discussed earlier, will cause the vector 6 to cone about the 
polar axis of the planet. The combined effect of these two variations gives rise to 
the change in P angle. 

High Inclination, Low-Earth Orbit 

This example illustrates a typical high inclination, low altitude Earth orbit. In this 
case, the analysis begins the first day of spring 1994. The complete list of orbital 
parameters used to initialize the problem are listed in Table 1. 
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Table 1. Orbital Parameters of Sample Case 1 


Semimajor Axis 
Eccentricity 

6785.58 km 
0.0 

Inclination 

51.6° 

Arg. of Pericenter 

Undefined 

Apsidal Rotation Rate 

3.7326° /day 

Initial Right Ascension 

358.77° 

Ascending Node Rate 

-4.99° /day 

Initial Solar Right Ascension 

0.5866° 

Initial Solar Declination 

0.2400° 


This case was selected to test the code over a beta angle range which would 
provide for the full range of shadowing situations. The orbit inclination and 
altitude selected provide for numerous precession cycles throughout the year as 
well as a number of periods of 100% sunlit orbits. 

The results of this case are shown in Figures 9 and 10 for the beta angle variation 
and time spent in shadow as a function of time. 


Shadow Terminator Calculations 
Sample Case 1 
(Beta Angle vs. Time) 



Figure 9. Case 1: Beta angle variation with time. 
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Shadow Terminator Calculation 
Sample Case 1 

(Percent time in Umbral and Penumbral Shadow) 



Figure 10. Case 1: Percent time in shadow. 

Sun-Synchronous Orbit 

A typical sun-synchronous orbit is modeled in this example. Sun-synchronous 
orbits are useful for mapping spacecraft since they are flown at nearly a constant 
beta angle. This provides a consistent lighting environment for the sunlit side of 
the orbit. This is accomplished through use of a retrograde (i.e., i > 90°) orbit 
which causes the ascending node to precess eastward. The altitude and inclination 
are selected such that the rate of movement of the ascending node closely 
matches the mean motion of the sun as it moves about the celestial sphere. A 
sun-synchronous orbit then should have a relatively flat beta angle versus time 
profile and an even flatter shadow time profile (since umbral shadow time 

variation is not extreme for a large range of beta angles about (5=0°). This example 
also tests the retrograde motion capability of the algorithm. 

The orbital parameters are presented in Table 2. The initial orbital parameters also 
correspond to the first day of spring 1994. The beta angle and shadow profiles are 
presented in Figures 11 and 12, respectively. 
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Beta Angle (deg.) 


Table 2. Orbital Parameters of Sample Case 2 


|| Semimajor Axis [ 

7083.14 km 


0.0 

Inclination 

98.2° 

Arg. of Pericenter 

Undefined 

Apsidal Rotation Rate 


Initial Right Ascension 

358.77° 1 

Ascending Node Rate 

| 

Initial Solar Right Ascension 

0.5866° 

Initial Solar Declination 

0.2400° 


Shadow Terminator Calculations 
Sample Case 2 
(Beta Angle vs. Time) 



Figure 11. Case 2: Beta angle variation versus time. 
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Figure 12. Case 2: Percent time in shadow. 


High Inclination Elliptical Orbit 

Sample case 3 tests the algorithm on a high inclination elliptical orbit. The 
inclination selected causes the apsidal rotation rate to go to zero. The benefits of 
this orbit inclination are exploited by the Molniya spacecraft which maintains the 
apogee and perigee at desired locations to facilitate communications. The elliptic 
nature of the orbit tests the robustness of the algorithm over a variety of 
conditions. Since the precession of the ascending node is slow (compared to a low 
altitude orbit) fewer cycles of beta are seen. 

The orbital parameters, corresponding to the first day of spring 1994, are presented 
in Table 3. The beta angle and shadow profiles are presented in Figures 13 and 14, 
respectively. 
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Beta Angle (deg.) 


Table 3. Orbital Parameters of Sample Case 3 


| Semimajor Axis 

42238.84 km 


0.8273 

Inclination 

63.3° 

Initial Arg. of Pericenter 

-63.3° 

Apsidal Rotation Rate 


Initial Right Ascension 

358.77° || 

Ascending Node Rate 

teEnssBWHBMI 

Initial Solar Right Ascension 

0.5866° 

Initial Solar Declination 

0.2400° 


Shadow Terminator Calculations 
Sample Case 3 
(Beta Angle vs. Time) 



Figure 13. Case 3: Beta angle variation versus time. 


23 















Shadow Terminator Calculation 
Sample Case 3 

(Percent time in Umbral and Penumbral Shadow) 



Figure 14. Case 3: Percent time in shadow. 


Geostationary Orbit 

A geostationary orbit is modeled in this sample case. The altitude is selected such 
that the orbit period matches the rotation rate of the planet. This has the effect of 
keeping the spacecraft over a given point on the planet. 

The orbital parameters are presented in Table 4. The beta angle and shadow 
profiles are presented in Figures 15 and 16, respectively. 


Table 4. Orbital Parameters of Sample Case 4 


Semimajor Axis 

42305.08 km 

Eccentricity 

0.00 

Inclination 

0.00° 

Arg. of Pericenter 

Undefined 

Apsidal Rotation Rate 

liWjKIMBBi 

Initial Right Ascension 

219.77° | 

Ascending Node Rate 

EEBIH 

Initial Solar Right Ascension 

0.5866° 

Initial Solar Declination 

0.2400° 
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Shadow Terminator Calculations 
Sample Case 4 
(Beta Angle vs. Time) 








CRRES Orbit 


This sample case was selected as a cross-check of the TCADS software, with the 
sample case presented in Mullins. 10 Mullins describes a method for solving the 
terminator problem through solution of a quartic polynomial. The shadow 
profile curve was successfully re-created using the algorithm. The beta angle 
profile was also checked against results from the Thermal Synthesizer System 
software. 5 

The orbital parameters are presented in Table 5. The beta angle and shadow 
profiles are presented in Figures 17 and 18, respectively. In this particular case, the 
Mullins assumption of constant parameters (i.e., variation parameters described 
by equations 49 and 50 are only updated at the beginning of an analysis orbit and 
held constant for the orbit) appears to be a valid assumption, since no significant 
differences in calculated percent time in shade are observed. 


Table 5. Orbital Parameters of Sample Case 5 


Semimajor Axis 

24450 km 

Eccentricity 

0.725 

Inclination 

00 

b 

o 

Initial Arg. of Pericenter 

180° 

Apsidal Rotation Rate 

0.7064°/day 

Initial Right Ascension 

68° 

Ascending Node Rate 

-0.3812°/day 

[ Initial Solar Right Ascension 

83.041° 

| Initial Solar Declination 

23.27° 
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Percent Time in Shadow 


Shadow Terminator Calculations 
Sample Case 5, CRRES 
(Beta Angle vs. Time) 



0 150 300 450 600 750 

Time (Days) 


Figure 17. Case 5: Beta angle variation versus time. 


Shadow Terminator Calculation 
Sample Case 5, CRRES 

(Percent time in Umbral and Penumbral Shadow) 



0 150 300 450 600 750 

Time (Days) 


Figure 18. Case 5: Percent time in shadow. 
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Conical versus Cylindrical Assumption 


This test case was created as an illustration of the potential for error of using a 
cylindrical shadow assumption. An elliptical orbit was specified using the 
Thermal Synthesizer System software 5 such that its path would barely skim the 
cylindrical umbral shadow. Next, the same orbit was modeled using the TCADS 
algorithm. The orbit parameters used are given in Table 6. 


Table 6. Orbital Parameters of Sample Case 6 


Semimajor Axis 

44859.14 km 

Eccentricity 

0.8408 

Inclination 

4.47° 

Initial Arg. of Pericenter 

90° 

Apsidal Rotation Rate 


Initial Right Ascension 

-90° i 

Ascending Node Rate 

ES59HH 

Initial Solar Right Ascension 

0.131° 

Initial Solar Declination 

0.054° 


The cylindrical assumption applied to the parameters presented above produced 
an umbral shadow time of 0.955% (—15 min.)(using the Thermal Synthesizer 
System software). A cylindrical shadow assumption does not provide for any 
penumbral shadow. When calculated using the routine created for TCADS, it was 
determined that the spacecraft would spend 0% of the orbit in the umbral shadow 
and 7.60% of the orbit period in the penumbral shadow: almost 2 hours in less 
than full sunlight conditions. The implication of this is that quickly reacting 
spacecraft components will be affected by this reduction in solar flux. Hence, a 
more accurate characterization of the umbral and penumbral shadows will lead to 
a more accurate thermal analysis with fewer required simplifying assumptions. A 
comparison of the cylindrical versus conical assumption is given in Table 7. 


Table 7. Comparison of Conical and Cylindrical Shadow Assumptions 


Component 

Conical Assumption 
(Minutes) 

Cylindrical Assumption 
(Minutes) 

| Umbra 1 

6.0 

15.1 i 

| Penumbra 

119.8 

Not calculated 1 
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